Exploratory Analysis for Russian MoTR Reading Data

shhh <- suppressPackageStartupMessages # It's a library, so shhh!

shhh(library( mgcv ))
shhh(library(dplyr))
shhh(library(ggplot2))
shhh(library(lme4))
shhh(library(tidymv))
shhh(library(gamlss))
shhh(library(gsubfn))
shhh(library(lmerTest))
shhh(library(tidyverse))
shhh(library(boot))
shhh(library(rsample))
shhh(library(plotrix))
shhh(library(ggrepel))
shhh(library(mgcv))
shhh(library(brms))
shhh(library(bayesplot))
shhh(library(tidyr))
shhh(library(car))
shhh(library(HDInterval))
shhh(library(gridExtra))
shhh(library(posterior))
shhh(library(readxl))
shhh(library(stringr))
shhh(library(loo))

shhh(library(coda))
shhh(library(cmdstanr))
shhh(library(rstan))
shhh(library(rstantools))

rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
rstan_options(auto_write = TRUE)
theme_set(theme_bw())
options(digits=4)
options(scipen=999)
set.seed(444)

Read in ET Data

file_list <- list.files("/Users/cui/Documents/uzh/PhD/Projects/Russian_Agreement/russian_gender/ref/Eyetracking/", pattern = "*.xlsx", full.names = TRUE)
et_raw <- file_list %>%
  lapply(read_excel) %>%
  bind_rows()

# View(et_raw)

Read in MoTR Data

# The path to the data
data_path <- "./data/"
data_names <- list.files(data_path)

# Read in the data from each participant and add to the data frame
motr_df <- data.frame()
for(name in data_names){
  subj <- gsub("reader_", "", gsub("_reading_measures.csv", "", name))
  temp_df <- read.csv(paste0(data_path, "/", name)) %>% mutate(subj_id = subj)
  motr_df <- rbind(motr_df, temp_df)
} 

motr_df <- motr_df %>% mutate(word_len = nchar(word),
                              word_length = scale(word_len)[,1]) %>% 
  group_by(subj_id, item_id) %>%
  arrange(subj_id, item_id) %>%
  mutate(word_len_pre1 = lag(word_length, n = 1),
         word_len_pre2 = lag(word_length, n = 2)) %>%
  ungroup()


View(motr_df)

# Clean the data
clean_df <- motr_df %>%
  # filter(subj_id != 171) %>%   # acc = 0.8
  filter(! list %in% c(98, 99)) %>% # filter practice and filler items
  mutate(skip = ifelse(total_duration==0, 1, 0),
         FPReg = ifelse(gaze_duration==0, NA, FPReg),
         FPFix = ifelse(gaze_duration==0, NA, FPFix)) %>%
  filter(skip == 0) %>%
  
  gather(measure, value, 18:26) %>%
  mutate(tgt_zero = if_else(measure %in% c("first_duration", "gaze_duration", "go_past_time", "right_bounded_rt", "total_duration") & value == 0, F, T)) %>%
  filter(tgt_zero != F) %>%
  dplyr::select(-tgt_zero, -cond_id, -skip, -word_len) %>%
  mutate(item_id = as.factor(item_id),
         subj_id = as.factor(subj_id)) %>%
  spread(measure, value) %>%
  gather(measure, value, c("first_duration", "gaze_duration", "go_past_time", "right_bounded_rt", "total_duration")) %>%
  mutate(outlier = value > (mean(value, na.rm = TRUE) + 3 * sd(value, na.rm = TRUE))) %>%
  filter(outlier == FALSE) %>%
  dplyr::select(-outlier) %>%
  spread(measure, value) %>%
  gather(measure, value, 21:29) %>%

  mutate(cond = case_when(
    target_gender == "M" & gender_match == "Mis" & type == "stim_adj" ~ "a",
    target_gender == "M" & gender_match == "Mis" & type == "stim_verb" ~ "b",
     target_gender == "M" & gender_match == "Mis" & type == "stim_pred_adj" ~ "c",
    target_gender == "M" & gender_match == "Match" & type == "stim_adj" ~ "d",
    target_gender == "M" & gender_match == "Match" & type == "stim_verb" ~ "e",
    target_gender == "M" & gender_match == "Match" & type == "stim_pred_adj" ~ "f",
    target_gender == "F" & gender_match == "Mis" & type == "stim_adj" ~ "g",
    target_gender == "F" & gender_match == "Mis" & type == "stim_verb" ~ "h",
    target_gender == "F" & gender_match == "Mis" & type == "stim_pred_adj" ~ "i",
    target_gender == "F" & gender_match == "Match" & type == "stim_adj" ~ "j",
    target_gender == "F" & gender_match == "Match" & type == "stim_verb" ~ "k",
    target_gender == "F" & gender_match == "Match" & type == "stim_pred_adj" ~ "l",
    TRUE ~ NA_character_ # This is the default case if none of the above conditions are met
  )) %>%
  dplyr::select(-list, -part, -type_id, -orig_item_number, -case, -animacy, -response_true, -response_chosen, -correctness) #%>%
  # drop_na()

clean_df <- clean_df %>%
  mutate(word = str_replace_all(word, "\\.", "")) %>%
  rowwise() %>%
  mutate(log_freq = ifelse(word %in% et_raw$IA_LABEL, 
                           et_raw$lg_frequency[match(word, et_raw$IA_LABEL)], 
                           NA_real_)) %>%
  ungroup()

View(clean_df)

correctness <- motr_df %>% dplyr::select(item_id, cond_id, subj_id, correctness) %>%
  filter(correctness != 99) %>% 
  distinct()

correctness_summary <- correctness %>%
  group_by(subj_id) %>%
  summarise(mean_correctness = mean(correctness),
            sd_correctness = sd(correctness),
            count = n())

View(correctness_summary)   # only subj_id 171 get acc = 0.8; others all > 0.88

# write.csv(correctness_summary, "./stats/correctness_summary.csv", row.names = FALSE)

RESEARCH QUESTIONS:

1. Are RTs different in gender-match versus gender-mismatch sentences? ==> main effect of grammaticality (gender match or not)

2. Are RTs different in Masculine versus Feminine sentence conditions? ==> main effect of gender of target word

3. Are RTs affected by sentence type (whether different lexical categories of the agreeing element will make the processing more difficult or not)? –> ADJ(adj + pre_adj) v.s. VERB ==> main effect of lexical type of sentences.

4. Are RTs affected by sentence type (whether agreeing element instantiates internal v.s. external agreement will make a difference in processing difficulty)? –> internal (modifying adjective) v.s. external (verb or predicative adjective) ==> main effect of syntax type of sentences.

5. Does the grammaticality effect within each lexical sentence type differ from each other? –> Whether the effect of grammaticality depends on the lexical type of the sentence (ADJ? VERB?) ==> interaction between grammaticality and lexical sentence type

6.Does the grammaticality effect within each syntax sentence type differ from each other? –> Whether the effect of grammaticality depends on the syntax type of the sentence (internal? external?) ==> interaction between grammaticality and syntax sentence type

7. Does the (possible) difference in the sensitivity to the grammaticality manipulation of Masculine versus Feminine conditions differ between lexical sentence types (ADJ v.s. VERB)? ==> 3-way interaction between grammaticality, gender and lexical sentence type

8. Does the (possible) difference in the sensitivity to the grammaticality manipulation of Masculine versus Feminine conditions differ between syntax sentence types (internal v.s. external)? ==> 3-way interaction between grammaticality, gender and syntax sentence type

contrast coding

# check conditions
clean_df$cond <- factor(clean_df$cond)
summary(clean_df$cond)

   a    b    c    d    e    f    g    h    i    j    k    l 
5922 6255 5805 6309 6498 6156 6129 6309 7623 5562 5616 6912 

clean_df <- clean_df %>% 
  mutate(
    #--------------------- main effects ---------------------
    Gram = ifelse(cond %in% c('a', 'b', 'c', 'g', 'h', 'i'), 1/6, -1/6), # Main effect grammaticality 
    Gen = ifelse(cond %in% c('a','b','c','d','e', 'f'), 1/6, -1/6), # Main effect gender
    TypL = ifelse(cond %in% c('a','c','d','f', 'g', 'i', 'j', 'l'), 1/8, -1/4), # Main effect of sentence type (ap vs v)
    TypS = ifelse(cond %in% c('a', 'd', 'g', 'j'), 1/4, -1/8), # Main effect of sentence type (a vs pv)
    
    #--------------------- 2 way interection ---------------------
    # Gram_x_TypL = ifelse(cond %in% c('a', 'c', 'i'), 1/8,
    #                  ifelse(cond %in% c('d', 'f', 'l'), -1/8,
    #                    ifelse(cond %in% c('e', 'k'), 1/4, -1/4))), # Grammaticality x type (ap v)
    
    # Gram_x_TypS = ifelse(cond %in% c('e', 'f', 'k', 'l'), 1/8,
    #                   ifelse(cond %in% c('b', 'c', 'h', 'i'), -1/8,
    #                          ifelse(cond %in% c('a', 'g'), 1/4, -1/4))), # Grammaticality x type (ap v)
    
    Gram_x_TypL = ifelse(cond %in% c('a', 'c', 'g', 'i', 'e', 'k'), 1/2, -1/2), # Grammaticality x type (ap v)
    Gram_x_TypS = ifelse(cond %in% c('a', 'g', 'e', 'f', 'k', 'l'), 1/2, -1/2), # Grammaticality x type (a pv)

    Gram_TypL_M = ifelse(cond %in% c('a', 'c', 'e'), 1/2, 
                    ifelse(cond %in% c('b', 'd', 'f'), -1/2, 0)), # gram x typl(ap v)_M
    Gram_TypS_M = ifelse(cond %in% c('a', 'e', 'f'), 1/2, 
                ifelse(cond %in% c('b', 'c', 'd'), -1/2, 0)), # gram x typs(a pv)_M
    Gram_TypL_F = ifelse(cond %in% c('g', 'i', 'k'), 1/2, 
                ifelse(cond %in% c('h', 'j', 'l'), -1/2, 0)), # gram x typl(ap v)_F
    Gram_TypS_F = ifelse(cond %in% c('a', 'e', 'f'), 1/2, 
                ifelse(cond %in% c('b', 'c', 'd'), -1/2, 0)), # gram x typs(a pv)_F
    
    #--------------------- 3 way interection ---------------------
    Gram_x_Gen_x_TypL = ifelse(cond %in% c('a', 'c', 'e', 'h', 'j', 'l'), 1/2, -1/2), # gen x typ1(ap v) x gram
    Gram_x_Gen_x_TypS = ifelse(cond %in% c('a', 'e', 'f', 'h', 'i', 'j'), 1/2, -1/2), # gen x typ1(ap v) x gram
    
    #--------------------- Within grammaticality type effects ---------------------
    Typ_Mis = ifelse(cond %in% c('a', 'c', 'g', 'i'), 1/4,
              ifelse(cond %in% c('b', 'h'), -1/2, 0)),  # type_Mis
    Typ_Match = ifelse(cond %in% c('d', 'f', 'j', 'l'), 1/4,
                  ifelse(cond %in% c('e', 'k'), -1/2, 0))  # type_Match
  ) %>% spread(measure, value) %>%
  # filter(word_nr == 3)
  filter(AOI_id == "R3")
  
clean_df

clean_df$target_gender <- factor(clean_df$target_gender)    # F, M
clean_df$gender_match <- factor(clean_df$gender_match)      # Match, Mismatch
clean_df$type <- factor(clean_df$type)                      # adj, pre_adj, verb

X_C_bar1 <- contr.sum(2)
X_C_bar2 <- contr.sum(3)

contrasts(clean_df$target_gender)<- X_C_bar1
contrasts(clean_df$gender_match)<- X_C_bar1
contrasts(clean_df$type)<- X_C_bar2

## Check contrasts
# contrasts(clean_df$target_gender)
# contrasts(clean_df$gender_match)
# contrasts(clean_df$type)    # The contrast is for: h0 -> grand mean; h1 -> adj-grand mean; h2 -> pre-gran mean

fuchs_stats_df = data.frame()
measure_types = c("gaze_duration", "go_past_time", "total_duration", 
                  "FPReg", "RegIn_incl")

for (meas in measure_types){
  print(paste("Fitting model for:", meas))
  
  if (meas %in% c("gaze_duration", "go_past_time", "total_duration")){
      model <- clean_df %>% filter(!is.na(.data[[meas]]))  %>% 
        # mutate(log_meas = log(.data[[meas]])) %>%
        lmer(as.formula(paste("log(", meas, ") ~ gender_match * type + target_gender + word_length + word_len_pre1 + word_len_pre2 + log_freq + 
            (1 + gender_match | subj_id) + (1 | item_id)")), 
            data = ., REML = F)
      coefs <- summary(model)$coefficients
      temp_results <- data.frame(
        measure = meas,
        beta = c("b_0", "b_Gram", "b_adj", "b_pred_adj", "b_fem", "b_len", "b_len_pre1", "b_len_pre2", "b_freq", "b_Gram_x_adj", "b_Gram_x_pred_adj"),
        bval = coefs[, "Estimate"],
        pval = coefs[, "Pr(>|t|)"]
      )
  }else{
      model <- clean_df %>% filter(!is.na(.data[[meas]]))  %>% 
        glmer(as.formula(paste(meas, "~ gender_match * type +  (1 | subj_id)")), 
            data = ., family=binomial(link = "logit"))
      coefs <- summary(model)$coefficients
      temp_results <- data.frame(
        measure = meas,
        beta = c("b_0", "b_Gram", "b_adj", "b_pred_adj", "b_Gram_x_adj", "b_Gram_x_pred_adj"),
        bval = coefs[, "Estimate"],
        pval = coefs[, "Pr(>|z|)"]
        )
  }
    fuchs_stats_df = rbind(fuchs_stats_df, temp_results)
}

fuchs_stats_df = fuchs_stats_df %>%
  mutate(sig = if_else(pval < 0.05, "SIG", ifelse(pval < .1, ".", "")))

# View(fuchs_stats_df)
# write.csv(fuchs_stats_df, "./stats/fuchs_replicates.csv", row.names = FALSE)

stats_df = data.frame()
measure_types = c("first_duration", "gaze_duration", "go_past_time", "total_duration", 
                  "FPReg", "RegIn_incl"
                  )

for (meas in measure_types){
  print(paste("Fitting model for:", meas))
  
  if (meas %in% c("first_duration", "gaze_duration", "go_past_time", "total_duration")){
      model <- clean_df %>% filter(!is.na(.data[[meas]]))  %>% 
        lmer(as.formula(paste("log(", meas, ") ~ Gram + Gen + TypL + TypS + Gram_x_TypL + Gram_x_TypS + Gram_x_Gen_x_TypL + Gram_x_Gen_x_TypS + 
            (1 | item_id) + (1 + Gram | subj_id)")), 
            data = ., REML = F)
      coefs <- summary(model)$coefficients
      temp_results <- data.frame(
        measure = meas,
        beta = c("b_0", "b_Gram", "b_Gen", "b_TypL", "b_TypS", 
                 "b_Gram_x_TypL", "b_Gram_x_TypS", "b_Gram_x_Gen_x_TypL", "b_Gram_x_Gen_x_TypS"),
        bval = coefs[, "Estimate"],
        pval = coefs[, "Pr(>|t|)"]
      )
  }else{
      model <- clean_df %>% filter(!is.na(.data[[meas]]))  %>% 
        glmer(as.formula(paste(meas, "~ Gram + Gram_x_TypL + Gram_x_TypS + Gram_x_Gen_x_TypL + Gram_x_Gen_x_TypS + 
            (1 | item_id) + (1 | subj_id)")), 
            data = ., family=binomial(link = "logit"))
      coefs <- summary(model)$coefficients
      temp_results <- data.frame(
        measure = meas,
        beta = c("b_0", "b_Gram", 
                 "b_Gram_x_TypL", "b_Gram_x_TypS", "b_Gram_x_Gen_x_TypL", "b_Gram_x_Gen_x_TypS"),
        bval = coefs[, "Estimate"],
        pval = coefs[, "Pr(>|z|)"]
        )
  }
    stats_df = rbind(stats_df, temp_results)
}

stats_df = stats_df %>%
  mutate(sig = if_else(pval < 0.05, "SIG", ifelse(pval < .1, ".", "")))

View(stats_df)

# write.csv(stats_df, "./stats/stats_lmer_new.csv", row.names = FALSE)

Model comparison

model_comparison_df = data.frame()
measure_types = c("gaze_duration", "go_past_time", "total_duration", 
                  "FPReg", "RegIn_incl")

for (meas in measure_types){
  print(paste("Compare models for:", meas))
    if (meas %in% c("gaze_duration", "go_past_time", "total_duration")){
        model_l <- clean_df %>% filter(!is.na(.data[[meas]])) %>%
        lmer(as.formula(paste("log(", meas, ") ~  Gram + Gen + TypL  + Gram_x_TypL + Gram_x_Gen_x_TypL + 
            (1 | item_id) + (1 + Gram | subj_id)")), data = ., REML = F)

        model_s <- clean_df %>% filter(!is.na(.data[["go_past_time"]]))  %>% 
        lmer(as.formula(paste("log(", meas, ") ~ Gram + Gen + TypS  + Gram_x_TypS + Gram_x_Gen_x_TypS + 
            (1 | item_id) + (1 + Gram | subj_id)")), data = ., REML = F)
    }else{
        model_l <- clean_df %>% filter(!is.na(.data[[meas]]))  %>% 
        glmer(as.formula(paste(meas, "~ Gram + Gram_x_TypL + Gram_x_Gen_x_TypL +
            (1 | item_id)")), 
            data = ., family=binomial(link = "logit"))
        model_s <- clean_df %>% filter(!is.na(.data[[meas]]))  %>% 
        glmer(as.formula(paste(meas, "~ Gram + Gram_x_TypS + Gram_x_Gen_x_TypS +
            (1 | item_id)")), 
            data = ., family=binomial(link = "logit"))
    }
  
    aic_bic_comparison <- data.frame(
      `Dependent Variable` = meas,
      `Model Type` = c("L model", "S model"),
      AIC = c(AIC(model_l), AIC(model_s)),
      BIC = c(BIC(model_l), BIC(model_s)))
    model_comparison_df = rbind(model_comparison_df, aic_bic_comparison)
}
# write.csv(model_comparison_df, "./stats/model_comparison_lmer.csv", row.names = FALSE)

Fit Bayesian models

function for creating stan data format

createStanDat<-function(d, dv,form){
  
  subj <- as.integer(factor(d$subj_id))
  N_subj <- length(unique(subj))
  item <- as.integer(factor(d$item_id))
  N_items <- length(unique(item))
  X <- unname(model.matrix(form, d))  
  attr(X, which="assign") <- NULL
  
  stanDat <- list(N = nrow(X),           
                  P = ncol(X),              
                  n_u = ncol(X),             
                  n_w = ncol(X),            
                  X = X,                     
                  Z_u = X,                 
                  Z_w = X,                   
                  J = N_subj,                
                  K = N_items,
                  dv = dv,                    
                  subj = subj,
                  item = item)
  stanDat
}

stan_gd <- createStanDat(d=subset(clean_df, !is.na(clean_df$gaze_duration)),
                             form=as.formula("~1+Gram+Gen+TypL+TypS+Gram_x_TypL+Gram_x_TypS+Gram_x_Gen_x_TypL+Gram_x_Gen_x_TypS"), 
                             dv=subset(clean_df, !is.na(clean_df$gaze_duration))$gaze_duration)

m1_gd <- stan(file = "stan/maxModel1.stan", 
                data = stan_gd,
                iter = 4000, 
                chains = 4,
                control = list(adapt_delta=0.99))
saveRDS(m1_gd, file = "model/m1_gd2.rds")

stan_gpt <- createStanDat(d=subset(clean_df, !is.na(clean_df$go_past_time)),
                             form=as.formula("~1+Gram+Gen+TypL+TypS+Gram_x_TypL+Gram_x_TypS+Gram_x_Gen_x_TypL+Gram_x_Gen_x_TypS"), 
                             dv=subset(clean_df, !is.na(clean_df$go_past_time))$go_past_time)

# sample from posterior distribution.
m1_gpt <- stan(file = "stan/maxModel1.stan", 
                data = stan_gpt,
                iter = 4000, 
                chains = 4,
                control = list(adapt_delta=0.99))
saveRDS(m1_gpt, file = "model/m1_gpt2.rds")

stan_td <- createStanDat(d=subset(clean_df, !is.na(clean_df$total_duration)),
                             form=as.formula("~1+Gram+Gen+TypL+TypS+Gram_x_TypL+Gram_x_TypS+Gram_x_Gen_x_TypL+Gram_x_Gen_x_TypS"), 
                             dv=subset(clean_df, !is.na(clean_df$total_duration))$total_duration)

# sample from posterior distribution.
m1_td <- stan(file = "stan/maxModel1.stan", 
                data = stan_td,
                iter = 4000, 
                chains = 4,
                control = list(adapt_delta=0.99))
saveRDS(m1_td, file = "model/m1_td2.rds")

binary dv

stan_fpreg <- createStanDat(d=subset(clean_df, !is.na(clean_df$FPReg)),
                             form=as.formula("~1+Gram+Gen+TypL+TypS+Gram_x_TypL+Gram_x_TypS+Gram_x_Gen_x_TypL+Gram_x_Gen_x_TypS"), 
                             dv=subset(clean_df, !is.na(clean_df$FPReg))$FPReg)

# sample from posterior distribution.
m1_fpreg <- stan(file = "stan/logitModel1.stan", 
                data = stan_fpreg,
                iter = 4000, 
                chains = 4,
                control = list(adapt_delta=0.99))
saveRDS(m1_fpreg, file = "model/m1_fpreg2.rds")

stan_regin <- createStanDat(d=subset(clean_df, !is.na(clean_df$RegIn_incl)),
                             form=as.formula("~1+Gram+Gen+TypL+TypS+Gram_x_TypL+Gram_x_TypS+Gram_x_Gen_x_TypL+Gram_x_Gen_x_TypS"), 
                             dv=subset(clean_df, !is.na(clean_df$RegIn_incl))$RegIn_incl)

# sample from posterior distribution.
m1_regin <- stan(file = "stan/logitModel1.stan", 
                data = stan_regin,
                iter = 4000, 
                chains = 4,
                control = list(adapt_delta=0.99))
saveRDS(m1_regin, file = "model/m1_regin2.rds")

Examine fitted stan models

# change m1_gd.rds to other models to check them

m1_gd <- readRDS("model/m1_fpreg.rds")

# check params
summary(m1_gd, pars = c('beta[1]', 'beta[2]', 'beta[3]', 'beta[4]', 'beta[5]', 'beta[6]', 'beta[7]', 'beta[8]', 'beta[9]'))

# check convergence
traceplot(m1_gd, pars = c("beta"))

# check predicts --> posterior parameter distr.
y_posterior <- extract(m1_gd) 
y_posterior$beta[,1] #intercept
y_posterior$beta[,2] #Gram
y_posterior$beta[,3] #Gen
y_posterior$beta[,4] #TypL
y_posterior$beta[,5] #TypS
y_posterior$beta[,6] #Gram_x_TypL
y_posterior$beta[,7] #Gram_x_TypS

# check predicts --> posterior parameter distr. back in normal space
pst_gram <- y_posterior$Gram
pst_gram
density_gram <- density(pst_gram)
plot(density_plot, main = "Density Plot of pst_gram", xlab = "pst_gram values", ylab = "Density", col = "red")

pst_gen <- y_posterior$Gen
pst_gen
density_gen <- density(pst_gen)
plot(density_gen, main = "Density Plot of pst_gen", xlab = "pst_gen values", ylab = "Density", col = "red")

pst_typl <- y_posterior$TypL
pst_typl
density_typl <- density(pst_typl)
plot(density_typl, main = "Density Plot of pst_typl", xlab = "pst_typl values", ylab = "Density", col = "red")

pst_typs <- y_posterior$TypS
pst_typs
density_typs <- density(pst_typs)
plot(density_typs, main = "Density Plot of pst_typs", xlab = "pst_typs values", ylab = "Density", col = "red")

pst_gramxtypl <- y_posterior$Gram_x_TypL
pst_gramxtypl
density_gramxtypl <- density(pst_gramxtypl)
plot(density_gramxtypl, main = "Density Plot of pst_gramxtypl", xlab = "pst_gramxtypl values", ylab = "Density", col = "red")

pst_gramxtyps <- y_posterior$Gram_x_TypS
pst_gramxtyps
density_gramxtyps <- density(pst_gramxtyps)
plot(density_gramxtyps, main = "Density Plot of pst_gramxtyps", xlab = "pst_gramxtyps values", ylab = "Density", col = "red")

pst_gramxgenxtypl <- y_posterior$Gram_x_Gen_x_TypL
pst_gramxgenxtypl
density_gramxgenxtypl <- density(pst_gramxgenxtypl)
plot(density_gramxgenxtypl, main = "Density Plot of pst_gramxgenxtypl", xlab = "pst_gramxgenxtypl values", ylab = "Density", col = "red")


pst_gramxgenxtyps <- y_posterior$Gram_x_Gen_x_TypS
pst_gramxgenxtyps
density_gramxgenxtyps <- density(pst_gramxgenxtyps)
plot(density_gramxgenxtyps, main = "Density Plot of pst_gramxgenxtyps", xlab = "pst_gramxgenxtyps values", ylab = "Density", col = "red")


# check posterior predicts --> the fits looks very good --> suspision of overfit? --> cv validation?
predicts <- y_posterior$Predict_rt
dim(predicts) # 8000 x 1271
y_true <- subset(clean_df, !is.na(clean_df$FPReg))$FPReg
ppc_dens_overlay(y_true, yrep = predicts[1:200, ])

calc_hpdi <- function(samples, prob = 0.95) {
  if(length(samples) < 2) {  # Check if there are too few samples to calculate an interval
    return(c(lower = NA, upper = NA))
  }
  # Ensure samples can be converted to an mcmc object
  tryCatch({
    hpd_interval <- HPDinterval(as.mcmc(samples), prob = prob)
    return(c(lower = hpd_interval[1], upper = hpd_interval[2]))
  }, error = function(e) {
    return(c(lower = NA, upper = NA))  # Return NA if any error occurs
  })
}

compile model results

stats_df <- data.frame()

# Define the measure types corresponding to your models
measure_types <- c("gaze_duration", "go_past_time", "total_duration", "FPReg", "RegIn_incl")

# Loop over each measure type to read the corresponding model and extract data
for (meas in measure_types) {
  if (meas == "gaze_duration") {
      m1 <- readRDS("model/m1_gd2.rds")
    } else if (meas == "go_past_time") {
      m1 <- readRDS("model/m1_gpt2.rds")
    } else if (meas == "total_duration") {
      m1 <- readRDS("model/m1_td2.rds")
    } else if (meas == "FPReg") {
      m1 <- readRDS("model/m1_fpreg2.rds")
    } else if (meas == "RegIn_incl") {
      m1 <- readRDS("model/m1_regin2.rds")
    }

  # Extract posterior distributions
  y_posterior <- extract(m1)
  intercept <- exp(y_posterior$beta[,1])

  betas <- c("b_0", "b_Gram", "b_Gen", "b_TypL", "b_TypS", "b_Gram_x_TypL", "b_Gram_x_TypS", "b_Gram_x_Gen_x_TypL", "b_Gram_x_Gen_x_TypS")
  posterior_samples <- list(intercept, y_posterior$Gram, y_posterior$Gen, y_posterior$TypL, y_posterior$TypS, y_posterior$Gram_x_TypL, y_posterior$Gram_x_TypS, y_posterior$Gram_x_Gen_x_TypL, y_posterior$Gram_x_Gen_x_TypS)
  
  hpdi_95 <- lapply(posterior_samples, function(x) hdi(x, credMass = 0.95))
  hpdi_89 <- lapply(posterior_samples, function(x) hdi(x, credMass = 0.89))

  # Prepare the results data frame
  temp_results <- data.frame(
    measure = rep(meas, length(betas)),
    beta = betas,
    bval_mean = sapply(posterior_samples, mean),
    crI_95_lower = sapply(posterior_samples, function(x) quantile(x, 0.025)),
    crI_95_upper = sapply(posterior_samples, function(x) quantile(x, 0.975)),
    crl_89_lower = sapply(posterior_samples, function(x) quantile(x, 0.055)),
    crl_89_upper = sapply(posterior_samples, function(x) quantile(x, 0.945)),
    hpdi_95_lower = sapply(hpdi_95, function(x) x[1]),
    hpdi_95_upper = sapply(hpdi_95, function(x) x[2]),
    hpdi_89_lower = sapply(hpdi_89, function(x) x[1]),
    hpdi_89_upper = sapply(hpdi_89, function(x) x[2]),
    bval_median = sapply(posterior_samples, median)
  )

  # Append the temp_results to the stats_df data frame
  stats_df <- rbind(stats_df, temp_results)
}

write.csv(stats_df, "./stats/stats_bayesian_new.csv", row.names = FALSE)

OBSERVATIONS:

gpt, gd, td, reg:

1. Effects of Grammaticality are always significant. –> There are main effect of grammaticality (gender match or not).

2. Effects of condTypS are significant or nearly significant. –> Very likely, there are main effect of syntactic agreement type (adj vs verb & pred_adj)

3. The significance of condTypL is always smaller than condTypS. In Bayesian, the crI of beta for TypL (diff between adj vs verb) is always larger and 0 is more to the middle of its distribution.

• If do sum contrast for each level of type, pred_adj type is significantly different from the grand mean.

4. For gd, the interaction between TypS and Gram is significant. For td and gpt, it is also near to significance (or 0 to be in the narrow tail of its distr. in Bayesian). –> external or internal agreement will affect the process of mismatches in sentence.

5. Interaction between TypL and Gram is far from significance.

6. In td, there are three way interaction between TypL & Gram & Gen. –> coincidence?

CONCLUSION:

1. Longer rt and more regressions in error setences.
2. Longer rt and possibly more regressions in external agreement than internal ones. –> external is more difficult.
3. But rt difference between mismatch and match is bigger in internal than in external. –> internal agreement kind of amplify the processing difficulty in mismatches. Maybe because external sentences are already difficult, so when combined with errors, the rt is not that different from in correct sentence. Maybe because people have their upboundary for the times they spend on difficult sentences? –> obviously, rt(difficult_external_type + error_sentence) < rt(difficult_external_type_sentence) + rt(error_sentence)

Bayesian model comparison

# Exclude outlier points 
clean_data_lexical <- clean_df %>% filter(subj_id != 131)
clean_data_syntax <- clean_df %>% filter(subj_id != 131)

# Refit the models on the cleaned data
lexical_gpt_clean <- createStanDat(d=subset(clean_data_lexical, !is.na(clean_data_lexical$go_past_time)),
                             form=as.formula("~1+Gram+Gen+TypL+Gram_x_TypL+Gram_x_Gen_x_TypL"), 
                             dv=subset(clean_data_lexical, !is.na(clean_data_lexical$go_past_time))$go_past_time)

# lexical_gpt
# sample from posterior distribution.
lexical_gpt_clean <- stan(file = "stan/Modelcomparison.stan", 
                data = lexical_gpt_clean,
                iter = 4000, 
                chains = 4,
                control = list(adapt_delta=0.99)
                )
# saveRDS(lexical_gpt_clean, file = "model/model_comparison/lexical_gpt_clean.rds")

syntax_gpt_clean <- createStanDat(d=subset(clean_data_syntax, !is.na(clean_data_syntax$go_past_time)),
                             form=as.formula("~1+Gram+Gen+TypS+Gram_x_TypS+Gram_x_Gen_x_TypS"), 
                             dv=subset(clean_data_syntax, !is.na(clean_data_syntax$go_past_time))$go_past_time)

# lexical_gpt
# sample from posterior distribution.
syntax_gpt_clean <- stan(file = "stan/Modelcomparison.stan", 
                data = syntax_gpt_clean,
                iter = 4000, 
                chains = 4,
                control = list(adapt_delta=0.99))
# saveRDS(syntax_gpt_clean, file = "model/model_comparison/syntax_gpt_clean.rds")

# Perform LOO on cleaned models
loo_lexical_clean <- loo(lexical_gpt_clean)
loo_syntax_clean <- loo(syntax_gpt_clean)

# Compare the cleaned models
loo_compare(loo_lexical_clean, loo_syntax_clean)

PLOT

agg_all = df %>%
  group_by(type, gender_match, region, measure) %>%
    summarise(
      m = mean(value),
      s = std.error(value),
      lower = m - 1.96 * s,
      upper = m + 1.96 * s
    ) %>%
  ungroup()

`summarise()` has grouped output by 'type', 'gender_match', 'region'. You can override using the `.groups` argument.

agg_all

fuchs_stats_df <-  read.csv("./stats/fuchs_replicates.csv")

Warning: cannot open file './stats/fuchs_replicates.csv': No such file or directoryError in file(file, "rt") : cannot open the connection

# Plot the reading times in each region broken down by stimulus type, reading measure and target gender

agg_all %>%
  # filter(type_s %in% c("stim_adj", "stim_verb", "stim_pred_adj")) %>%
  filter(measure %in% c("gaze_duration", "go_past_time", "total_duration")) %>%
  
  mutate(type = factor(type, levels = c("stim_adj", "stim_pred_adj", "stim_verb"), labels=c("Modifying Adj.", "Predicative Adj.",  "Verb"))) %>%
  mutate(measure = factor(measure, levels = c("gaze_duration", "go_past_time", "total_duration"), labels=c("Gaze Duration", "Go Past Times", "Total Duration"))) %>%
  
  ggplot(aes(x=region, y=m, color = gender_match)) +
    geom_rect(aes(xmin=2.5, xmax=3.5, ymin=lower-100, ymax=upper+100), color = NA, fill = "green", alpha=0.01) +
    geom_point() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width = 0.2) +
    geom_line() +
    facet_grid(measure~type, scales = "free_y") +
    # scale_x_continuous(breaks = 1:4, labels = c("critical -2", "critical -1", "critical", "critical +1")) +
    #ggtitle(paste0("Region by region plot for ", stim)) +
    ylab("Reading time (ms)") +
    xlab("Sentence Region") +
  scale_color_manual(values = c("Match" = "#56BCC2", "Mis" = "#E77D72")) +
  #coord_cartesian(ylim=c(0, 900)) +
  theme(
    # axis.text.x = element_text(size = 9),
    legend.position = "bottom"
  )

ggsave(paste0("./images/RT_results_all_types.pdf"), device="pdf", height=5, width=5)

# Plot the reading times in each region broken down by stimulus type, reading measure and target gender

agg_s %>%
  filter(measure %in% c("FPReg", "RegIn_incl")) %>%
  
  mutate(type_s = factor(type_s, levels = c("Modifying_Adj", "Predictive_Adj_Verb"), labels=c("Modifying Adj.", "Predicative Adj. & Verb"))) %>%
    mutate(measure = factor(measure, levels = c("FPReg", "RegIn_incl"), labels=c("Regr. out Prob.", "Regr. in Prob."))) %>%

  ggplot(aes(x=region, y=m, color = gender_match)) +
    geom_rect(aes(xmin=2.5, xmax=3.5, ymin=lower-0.2, ymax=upper+0.2), color = NA, fill = "green", alpha=0.01) +
    geom_point() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width = 0.2) +
    geom_line() +
    facet_grid(measure~type_s, scales = "free_y") +
    ylab("Regression Prob.") +
    xlab("Sentence Region") +
  scale_color_manual(values = c("Match" = "#56BCC2", "Mis" = "#E77D72")) +
  #coord_cartesian(ylim=c(0, 900)) +
  theme(
    # axis.text.x = element_text(size = 9),
    legend.position = "bottom"
  )
ggsave(paste0("./images/Regression_results_synt_position2.pdf"), device="pdf", height=3.4, width=5)

agg_s2 = df %>%
  group_by(type_s, subj_id, region, measure) %>%
    summarise(
      m = mean(value),
    ) %>%
  ungroup() %>%
  group_by(type_s, region, measure) %>%
  pivot_wider(
    names_from = type_s,
    values_from = m,
    names_prefix = "mean_"
  ) %>%
  # Calculate the difference between 'Mis' and 'Match'
  drop_na() %>%
  mutate(
    diff = mean_Predictive_Adj_Verb - mean_Modifying_Adj
  ) %>%
  group_by(region, measure) %>%
  summarise(
    m_diff = mean(diff),
    s = std.error(diff),
    lower = m_diff - 1.96 * s,
    upper = m_diff + 1.96 * s
  ) %>%
  ungroup()

`summarise()` has grouped output by 'type_s', 'subj_id', 'region'. You can override using the `.groups` argument.`summarise()` has grouped output by 'region'. You can override using the `.groups` argument.

agg_s3 = df %>%
  group_by(type_s, gender_match, item_id, region, measure) %>%
    summarise(
      m = mean(value)
    ) %>%
  ungroup() %>%
  group_by(type_s, region, measure) %>%
  pivot_wider(
    names_from = gender_match,
    values_from = m,
    names_prefix = "mean_"
  ) %>%
  # Calculate the difference between 'Mis' and 'Match'
  drop_na() %>%
  mutate(
    diff = mean_Mis - mean_Match
  ) %>%
  group_by(type_s, region, measure) %>%
  summarise(
    m_diff = mean(diff),
    s = std.error(diff),
    lower = m_diff - 1.96 * s,
    upper = m_diff + 1.96 * s
  ) %>%
  ungroup()

`summarise()` has grouped output by 'type_s', 'gender_match', 'item_id', 'region'. You can override using the `.groups` argument.`summarise()` has grouped output by 'type_s', 'region'. You can override using the `.groups` argument.

agg_s3 %>%
  # filter(type_s %in% c("stim_adj", "stim_verb", "stim_pred_adj")) %>%
  filter(measure %in% c("FPReg", "RegIn_incl")) %>%
  
  mutate(type_s = factor(type_s, levels = c("Modifying_Adj", "Predictive_Adj_Verb"), labels=c("Modifying Adj.", "Predicative Adj. & Verb"))) %>%
  mutate(measure = factor(measure, levels = c("FPReg", "RegIn_incl"), labels=c("First Pass Regression out Prob.", "Regression in Prob."))) %>%

  ggplot(aes(x=region, y=m_diff, color = type_s)) +
    geom_rect(aes(xmin=2.5, xmax=3.5, ymin=-0.3, ymax=0.3), color = NA, fill = "green", alpha=0.01) +
    geom_point() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width = 0.2) +
    geom_line() +
    facet_wrap(~ measure, scales = "free_y") +
    #ggtitle(paste0("Region by region plot for ", stim)) +
    labs(
      title = "Regression Difference under Mismatch and Match conditions",
      y = "Regression prob. difference",
      x = "Sentence Region",
      color = "Agreement Type"
    ) +
  scale_x_continuous(breaks = c(1:5)) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5)
  )

ggsave(paste0("./images/Regression_results_interaction_in_synt.pdf"), device="pdf", height=5, width=8)

agg_l2 = df %>%
  group_by(type_l, subj_id, region, measure) %>%
    summarise(
      m = mean(value),
    ) %>%
  ungroup() %>%
  group_by(type_l, region, measure) %>%
  pivot_wider(
    names_from = type_l,
    values_from = m,
    names_prefix = "mean_"
  ) %>%
  # Calculate the difference between 'Mis' and 'Match'
  drop_na() %>%
  mutate(
    diff = mean_Modifying_Adj_Predictive_Adj - mean_Verb
  ) %>%
  group_by(region, measure) %>%
  summarise(
    m_diff = mean(diff),
    s = std.error(diff),
    lower = m_diff - 1.96 * s,
    upper = m_diff + 1.96 * s
  ) %>%
  ungroup()

`summarise()` has grouped output by 'type_l', 'subj_id', 'region'. You can override using the `.groups` argument.`summarise()` has grouped output by 'region'. You can override using the `.groups` argument.

# Plot the reading times in each region broken down by stimulus type, reading measure and target gender

agg_l2 %>%
  # filter(type_s %in% c("stim_adj", "stim_verb", "stim_pred_adj")) %>%
  filter(measure %in% c("gaze_duration", "go_past_time", "total_duration")) %>%
  mutate(measure = factor(measure, levels = c("gaze_duration", "go_past_time", "total_duration"), labels=c("Gaze Duration", "Go Past Times", "Total Duration"))
         ) %>%

  ggplot(aes(x=region, y=m_diff, color = measure)) +
    geom_rect(aes(xmin=2.5, xmax=3.5, ymin=lower-100, ymax=upper+100), color = NA, fill = "green", alpha=0.01) +
    geom_point() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width = 0.2) +
    geom_line() +
    facet_wrap(~ measure, scales = "free_y") +
    labs(
      title = "Difference in RTs:\n Adj. (Modifying Adj. and Predictive Adj.) vs Verb",
      y = "Reading time difference (ms)",
      x = "Sentence Region",
      color = "Measure"
    ) +
  #coord_cartesian(ylim=c(0, 900)) +
  scale_x_continuous(breaks = c(1:5)) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5)
  )

ggsave(paste0("./images/RT_results_main_diff_in_lexical.pdf"), device="pdf", height=5, width=8)

# Plot the reading times in each region broken down by stimulus type, reading measure and target gender

agg_l2 %>%
  # filter(type_s %in% c("stim_adj", "stim_verb", "stim_pred_adj")) %>%
  filter(measure %in% c("FPReg", "RegIn_incl")) %>%
  mutate(measure = factor(measure, levels = c("FPReg", "RegIn_incl"), labels=c("First Pass Regression out Prob.", "Regression in Prob."))) %>%

  ggplot(aes(x=region, y=m_diff, color = measure)) +
    geom_rect(aes(xmin=2.5, xmax=3.5, ymin=-0.2, ymax=0.2), color = NA, fill = "green", alpha=0.01) +
    geom_point() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width = 0.2) +
    geom_line() +
    facet_wrap(~ measure, scales = "free_y") +
    labs(
      title = "Difference in Regressions: \nAdj. (Modifying Adj. and Predictive Adj.) vs Verb",
      y = "Regression prob. difference",
      x = "Sentence Region",
      color = "Measure"
    ) +
  #coord_cartesian(ylim=c(0, 900)) +
  scale_x_continuous(breaks = c(1:5)) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5)
  )

ggsave(paste0("./images/Regression_results_main_diff_in_lexical.pdf"), device="pdf", height=5, width=8)

agg_l3 = df %>%
  group_by(type_l, gender_match, item_id, region, measure) %>%
    summarise(
      m = mean(value)
    ) %>%
  ungroup() %>%
  group_by(type_l, region, measure) %>%
  pivot_wider(
    names_from = gender_match,
    values_from = m,
    names_prefix = "mean_"
  ) %>%
  # Calculate the difference between 'Mis' and 'Match'
  drop_na() %>%
  mutate(
    diff = mean_Mis - mean_Match
  ) %>%
  group_by(type_l, region, measure) %>%
  summarise(
    m_diff = mean(diff),
    s = std.error(diff),
    lower = m_diff - 1.96 * s,
    upper = m_diff + 1.96 * s
  ) %>%
  ungroup()

`summarise()` has grouped output by 'type_l', 'gender_match', 'item_id', 'region'. You can override using the `.groups` argument.`summarise()` has grouped output by 'type_l', 'region'. You can override using the `.groups` argument.

# Plot the reading times in each region broken down by stimulus type, reading measure and target gender

agg_l3 %>%
  # filter(type_s %in% c("stim_adj", "stim_verb", "stim_pred_adj")) %>%
  filter(measure %in% c("gaze_duration", "go_past_time", "total_duration")) %>%
  
  mutate(type_l = factor(type_l, levels = c("Modifying_Adj_Predictive_Adj", "Verb"), labels=c("Modifying Adj. & Predicative Adj.", "Verb"))) %>%
  mutate(measure = factor(measure, levels = c("gaze_duration", "go_past_time", "total_duration"), labels=c("Gaze Duration", "Go Past Times", "Total Duration"))) %>%

  ggplot(aes(x=region, y=m_diff, color = type_l)) +
    geom_rect(aes(xmin=2.5, xmax=3.5, ymin=lower-100, ymax=upper+100), color = NA, fill = "green", alpha=0.01) +
    geom_point() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width = 0.2) +
    geom_line() +
    facet_wrap(~ measure, scales = "free_y") +
    #ggtitle(paste0("Region by region plot for ", stim)) +
    labs(
      title = "RTs Difference under Mismatch and Match conditions",
      y = "Reading time difference (ms)",
      x = "Sentence Region",
      color = "Agreement Type"
    ) +
  scale_x_continuous(breaks = c(1:5)) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5)
  )

ggsave(paste0("./images/RT_results_interaction_in_lexical.pdf"), device="pdf", height=5, width=8)

agg_l3 %>%
  # filter(type_s %in% c("stim_adj", "stim_verb", "stim_pred_adj")) %>%
  filter(measure %in% c("FPReg", "RegIn_incl")) %>%
  
  mutate(type_l = factor(type_l, levels = c("Modifying_Adj_Predictive_Adj", "Verb"), labels=c("Modifying Adj. & Predicative Adj.", "Verb"))) %>%
  mutate(measure = factor(measure, levels = c("FPReg", "RegIn_incl"), labels=c("First Pass Regression out Prob.", "Regression in Prob."))) %>%

  ggplot(aes(x=region, y=m_diff, color = type_l)) +
    geom_rect(aes(xmin=2.5, xmax=3.5, ymin=-0.3, ymax=0.3), color = NA, fill = "green", alpha=0.01) +
    geom_point() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width = 0.2) +
    geom_line() +
    facet_wrap(~ measure, scales = "free_y") +
    #ggtitle(paste0("Region by region plot for ", stim)) +
    labs(
      title = "Regression Difference under Mismatch and Match conditions",
      y = "Regression prob. difference",
      x = "Sentence Region",
      color = "Agreement Type"
    ) +
  scale_x_continuous(breaks = c(1:5)) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5)
  )

ggsave(paste0("./images/Regression_results_interaction_in_lexical.pdf"), device="pdf", height=5, width=8)

# Plot the reading times in each region broken down by stimulus type, reading measure and target gender

agg_sl2 %>%
  filter(measure %in% c("First Pass Regression out Prob.", "Regression in Prob.")) %>%
  ggplot(aes(x=region, y=m_diff, color = measure, group = Prediction, linetype = Prediction)) +
    geom_rect(aes(xmin=2.5, xmax=3.5, ymin=-0.2, ymax=0.2), color = NA, fill = "green", alpha=0.01) +
    geom_hline(yintercept = 0, color = "gray30") + 
    geom_point(aes(shape = Prediction), position = position_dodge(width = 0.3)) +
    geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, position = position_dodge(width = 0.3)) +
    geom_line(position = position_dodge(width = 0.3)) +
    facet_wrap(~ measure, scales = "free_y") +
    labs(
      title  = "Feature-matching Mechanism vs Lexical category:\nReg_Prob.(Agreement - Concord) vs Reg_Prob.(Adj - Verb)",
      y = "Regression difference (ms)",
      x = "Sentence Region",
      linetype = "Prediction"
    ) +
  #coord_cartesian(ylim=c(0, 900)) +
  scale_x_continuous(breaks = c(1:5)) +
  scale_linetype_manual(values = c("solid", "dashed")) +
  guides(color = FALSE, shape = FALSE) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5)
  )

ggsave(paste0("./images/Regression_main_diff_both_thoery.pdf"), device="pdf", height=4, width=16/3)

agg_sl3 %>%
  filter(measure %in% c("Gaze Duration", "Go Past Times", "Total Duration")) %>%
  ggplot(aes(x = region, y = m_diff, color = type, group = interaction(Prediction, type), linetype = Prediction)) +
    geom_rect(aes(xmin = 2.5, xmax = 3.5, ymin = lower - 100, ymax = upper + 100), color = NA, fill = "green", alpha = 0.01) +
    geom_hline(yintercept = 0, color = "gray30") + 
    geom_point(aes(shape = type)) +
    geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
    geom_line() +
    geom_text(aes(label = annotation, y = annotation_y), vjust = 0, color = "black") +
    facet_grid(Prediction ~ measure, scales = "free_y") +
    labs(
      title = "Interaction between Grammaticality and \n Feature-match Mechanism / Lexical Category",
      y = "Reading time difference (Mis - Match)",
      x = "Sentence Region"
    ) +
  scale_x_continuous(breaks = c(1:5)) +
  scale_color_manual(values = c(
    "Concord" = "#56BCC2",
    "Agreement" = "#E77D72",
    "Adjective" = "#56BCC2",
    "Verb" = "#E77D72"
  )) +
  scale_shape_manual(values = c(
    "Concord" = 16, # Filled circle
    "Agreement" = 17, # Filled triangle
    "Adjective" = 16, # Filled circle
    "Verb" = 17 # Filled triangle
  )) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5)
  ) +
  guides(
    linetype = "none",
    color = guide_legend(
      title = "Type",
      ncol = 4,
      byrow = TRUE,
      override.aes = list(
        linetype = c("solid", "solid", "dotdash", "dotdash"),
        shape = c(16, 17, 16, 17)
      )
    ),
    shape = guide_legend(
      title = "Type",
      ncol = 4,
      byrow = TRUE,
      override.aes = list(
        linetype = c("solid", "solid", "dotdash", "dotdash")
      )
    )
  )
ggsave(paste0("./images/RT_interaction_both_thoery.pdf"), device="pdf", height=5, width=8)

agg_sl3 %>%
  filter(measure %in% c("First Pass Regression out Prob.", "Regression in Prob.")) %>%
  ggplot(aes(x = region, y = m_diff, color = type, group = interaction(Prediction, type), linetype = Prediction)) +
    geom_rect(aes(xmin = 2.5, xmax = 3.5, ymin = lower - 0.2, ymax = upper + 0.2), color = NA, fill = "green", alpha = 0.01) +
    geom_hline(yintercept = 0, color = "gray30") + 
    geom_point(aes(shape = type)) +
    geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
    geom_line() +
    # geom_text(aes(label = annotation, y = annotation_y), vjust = 0, color = "black") +
    facet_grid(Prediction ~ measure, scales = "free_y") +
    labs(
      title = "Interaction between Grammaticality and \n Feature-match Mechanism / Lexical Category",
      y = "Regression prob. difference (Mis - Match)",
      x = "Sentence Region"
    ) +
  scale_x_continuous(breaks = c(1:5)) +
  scale_color_manual(values = c(
    "Concord" = "#56BCC2",
    "Agreement" = "#E77D72",
    "Adjective" = "#56BCC2",
    "Verb" = "#E77D72"
  )) +
  scale_shape_manual(values = c(
    "Concord" = 16, # Filled circle
    "Agreement" = 17, # Filled triangle
    "Adjective" = 16, # Filled circle
    "Verb" = 17 # Filled triangle
  )) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5)
  ) +
  guides(
    linetype = "none",
    color = guide_legend(
      title = "Type",
      ncol = 4,
      byrow = TRUE,
      override.aes = list(
        linetype = c("solid", "solid", "dotdash", "dotdash"),
        shape = c(16, 17, 16, 17)
      )
    ),
    shape = guide_legend(
      title = "Type",
      ncol = 4,
      byrow = TRUE,
      override.aes = list(
        linetype = c("solid", "solid", "dotdash", "dotdash")
      )
    )
  )
ggsave(paste0("./images/Regression_interaction_both_thoery.pdf"), device="pdf", height=5, width=16/3)